Exact Solution of a Reaction-Diffusion Model with Particle Number Conservation 
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We analytically investigate a Id branching-coalescing model with reflecting boundaries in a canon- 
ical ensemble where the total number of particles on the chain is conserved. Exact analytical calcu- 
lations show that the model has two different phases which are separated by a second-order phase 
transition. The thermodynamic behavior of the canonical partition function of the model has been 
calculated exactly in each phase. Density profiles of particles have also been obtained explicitly. 
It is shown that the exponential part of the density profiles decay on three different length scales 
which depend on total density of particles. 
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I. INTRODUCTION 



Recently much attention has been payed to the study 
of shocks in one-dimensional reaction-diffusion models 
[l|, 13, y> 13 • The shocks are defined as discontinuities in 
the space dependence of density of particles in the sys- 
tem and behave as collective excitations in system. They 
can be characterized by their position which performs 
a random walk. The best known example in which the 
shock can appear is asymmetric simple exclusion process 
(ASEP) with open boundaries Q. The mathematical 
relevance of the ASEP is that it is a discrete version of 
the Burgers equation in an appropriate scaling limit. The 
ASEP contains one class of particles (first class particles) 
which can be injected and extracted from the boundaries 
of a one-dimensional chain while hopping in the bulk with 
asymmetric rates. The ASEP has several applications to 
the realistic systems. For instance, it can be considered 
as a simple model for traffic flow lq|. 
There are different ways to provoke a shock in one- 
dimensional reaction-diffusion models. One can consider 
the ASEP on a closed chain in the presence of a sec- 
ond class particle. Compared to the first class parti- 
cles, the second class particles move very slowly. In 
nil, H E m in the shape of the shock is calculated 
as seen from a second class particle. Another method 
is to introduce a slow link in the system [13|. The first 
class particles cross this link with a smaller crossing rate 
than that of the other links in the system. In this case 
the width of the shock as a function of the length of the 
system L scales as L^/^ or L^/^ depending on whether 
the density of particles is equal to ^ or not [TJ • Shocks 
have also been observed in the ASEP with creation and 
annihilation of particles in the bulk of the system |15llla |. 
In a recent paper we have numerically studied the shocks 
in a spatially asymmetric one-dimensional branching- 
coalescing model with reflecting boundaries in a canon- 
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ical ensemble [I3. In this model the particles diffuse, 
coagulate and decoagulate on a lattice of length L; how- 
ever, the total number of the particles is kept fixed. It is 
predicted that the model has two different phases and in 
one phase the density profile of the particles has a shock 
structure. We have confirmed our numerical results by 
using the Yang-Lee theory of phase transitions [I^] which 
has recently been shown to be applicable to the study of 
critical behaviors of out-of-equilibrium systems 19, 2(f\. 
In the present work we will show that by working in the 
canonical ensemble, the model is exactly solvable in the 
sense that the thermodynamic limit of physical quanti- 
ties can be calculated exactly. The canonical partition 
function of the model defined as the sum over stationary 
state weights can also be calculated exactly. By applying 
the Yang-Lee theory we can calculate the line of the par- 
tition function zeros; and therefore, spot the transition 
point. The order of the transition can also be identi- 
fied by investigating the density of these zeros near the 
critical point. We will also obtain the exact expressions 
for the density profile of the particles on the chain in the 
thermodynamic limit. This paper is organized as follows: 
In Section II we will define the model and introduce the 
mathematical preliminaries. In Section III we will calcu- 
late the canonical partition function of our model using 
a matrix product formalism and find its behavior in the 
thermodynamic limit. We will also find the line of canon- 
ical partition function zeros to confirm our numerical re- 
sults in [131 • In Section IV we will calculate the density 
profile of the particles on the chain in each phase. In the 
last section we will discuss the results and compare them 
with the case where the total number of particles is not 
conserved. 



II. THE MODEL: MATHEMATICAL 
PRELIMINARIES 

In this section we will briefly review the definition of 
the model and also define its grand canonical partition 
function. We will then calculate the canonical partition 



function of the model explicitly. The model consists of 
one class of particles which diffuse on a one-dimensional 
chain of length L. Whenever two of these particles meet, 
they can coagulate to a single particle. In the same way, a 
single particle can decoagulate into two particles. There 
is no particle input or output at the boundaries. The 
reaction rules between two consecutive sites i and i + 1 
on the chain are explicitly as follows: 



(1) 



in which A and stand for the presence of a particle 
and a hole respectively. As can be seen, the parameter 
q determines the asymmetry of the model. For q > 1 
(q < 1) the particles have a tendency to move in the 
leftward (rightward) direction. For any q the model is 
also invariant under the following transformations: 
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in which i is a given site on the chain. One should also 
note that the rules (Q do not conserve the number of 
particles and therefore the model should be studied in a 
grand canonical ensemble. The model without particle 
number conservation has already been studied both us- 
ing Empty Interval Method (EIM) and Matrix Product 
Formahsm (MPF) [Hlg^ . It turns out that the model 
has two different phases in this case: Two exponential 
phases which are called the hight-density and the low- 
density phases. On the coexistence line of these phases 
the density of the particles on the chain has a linear pro- 
file. It is known that the phase diagram of the ASEP 
contains a first-order phase transition line where the in- 
jection and the extraction rates are equal and less than 
i . Along this line the density profile of particles is linear 
which is a consequence of superposition of states where a 
shock between a low-density region and a hight-density 
region is present at an arbitrary position |5|, l2l| . As we 
will see the linear profile in present model can also be 
interpreted as a sign for the existence of shocks. 
In order to study the shocks we restrict the total num- 
ber of particles on the chain to be M so that their to- 
tal density is always equal to p = ■^. This means that 
we are working in a canonical ensemble. The station- 
ary state probability distribution function can be cal- 
culated using the MPF as follows: We assign two non- 
commuting operators D and _E to a particle and a hole 
respectively. Now the probability for occurring any con- 
figuration {t} — {ti,--- ,ri} in the steady state with 
exactly M particles can be obtained from 



P{{r}) 



3(M-j:l,n) 



Zl. 



M 



{W\\{{nD + {l-n)E)\V) 

(3) 



in which r^ = 1 if the site i is occupied by a particles and 
Ti = if it is empty. The normalization factor Z^^m in 
©, which will be called the canonical partition function 
of the model hereafter, should be obtained from the fact 
that X]{t> P{{'''}) — 1- It '^^'^ be writen as 

L L 

Zlm = J2 ^(^'^ - J2 ^*)(^l n(^''^ + (1 " ri)E)\V). 

{t} i=l i=l 

(4) 
The Dirac delta in Q and Q guaranties the total num- 
ber of particles to be M in the steady state. In order to 
have a stationary probability distribution, the operators 
D and E besides the vectors \V) and {W\ should satisfy 
the following quadratic algebra [23 



[E, E] = q 

ED - ED = q{l 
DE -DE 



A)ED- \DE - iL>2 
-qED + ^DE - qD^ 
DD~ DD = -qAED ~ ^DE +{q+ ^)D^ 
{W\E = {W\D = , E\V) ^ D\V) ^ 0. 



(5) 



The operators D and E are auxiliary operators and do 
not enter into calculating ^ and Q . The following four- 
dimensional representation has been found for the alge- 
bra ® [H 



D 



E = 




(6) 



in which a and b are arbitrary constants and | W) is sim- 
ply transpose of {W\. The matrix representations for D 
and E are also given in [23 ■ Using ^ one can calculate 
the steady state weight of any given configuration. 
It turns out that the direct calculation of Q) is not always 
an easy task; therefore, we define the grand canonical 
partition function which can easily be calculated: 

L 
{r} ^=l 

= Y,e'z,,M (7) 

in which ^ is the fugacity associated with the particles. 
The total density of particles p should then be fixed by 
the fugacity of them through the following equation 



P 



lim 



Ld£, 



lnZi(0. 



(8) 



One can expect that each value of the fugacity ^ corre- 
sponds to each value of the total density. In this case, the 



density-fugacity relation ^ is invertible and the equiv- 
alence of the canonical and grand canonical ensemble 
holds. After calculating the grand canonical partition 
function iQ, one can invert the series to calculate the 
canonical partition function using 



Z 



L.M 



1 

27ri 



d£. 






(9) 



where C is a contour which encircles the origin anti- 
clockwise. For our model; however, there appears a situ- 



ation where the equivalence of ensembles fails in a special 
region in the parameters space. There is the place where 
the shocks appear in the system. 

As an important physical quantity one can study the den- 
sity profile of particles on the chain in the canonical en- 
semble; nevertheless, the calculation of the density pro- 
file of the particles is much more easily done in the grand 
canonical ensemble. Let us define the unnormalized av- 
erage particle number at site i in the grand canonical 
ensemble as: 



J 



{r} J = l 3=i+l 



(10) 



r 



We will then translate the results in the grand canonical 
ensemble into those in the canonical ensemble using the 
following formula 



\Pi/L.M 






(11) 



As in lO the contour C in (|ll|l encircles the origin anti- 
clockwise. In (|1(J|) and (|ll|l the superscript (u) means 
that it is an unnormalized quantity. The normalized av- 
erage particle number at site i should be obtained from 

{Pi) = {Pi)''lm/Zl,m- 



@ one can obtain the canonical partition function by ap- 
plying the steepest decent method. The grand canonical 
partition function of this model can easily be calculated 
using 10 and is simply given by: 



Zl{0 = {W\{iD + E)^\V) = {W\C'^\V) 



(12) 



in which we have defined C :~ ^D + E. The matrix 
representations for the operators D and E are given by 
©. After some algebra we find 



III. CANONICAL AND GRAND CANONICAL 
PARTITION FUNCTIONS 



.(1)/ 



zdo - zi^'io + zfio + ^r(o + ^r(o (is) 



.(3), 



7{i)i 



In this section we will calculate the grand canonical 
partition function of the model explicitly and then using 



in which 



J 



Z^l\0 = 



Z^l\0 



4''(0 



-q^^e 



](i + eA)^ 



^(<z2-(i + ^A))(g2(i + eA)-l 

r qH<f-l){l + ^^) 2L 

\q^ + l)(g2 - (1 + ^A))(g2(l + A) - (1 + ^A)) ^'^ 

. -g^(g^-l)(l + gA) 

\q^ + 1)((1 + A) - ^2(1 + ^A))(l - ^2(1 + ^A))^ 



'2L 



q^A{^ - 1) 



:]( 



l + CA,^ 



^(g2(i + A) - (1 + ^A)){q^l -h CA) - (l-f A))^^ 1 + A 

I 



(14) 
(15) 
(16) 
(17) 



Because of the symmetry of the model (0) one will only will consider the case g > 1 hereafter, and the results 
need to study either the case q > 1 or g < 1. We for the case q < 1 can easily be obtained by applying 



the transformations (0). Obviously for g > 1 we have 
q^ > q~^. On the other hand since A,^ > we always 
have (1 + ^A) > ( j^^ )■ Now two different cases can be 
distinguished: We will either have 1 < q < vT+"^S or 
1 < VI + ?^ < Q- For these two cases the asymptotic 
behaviors of the grand canonical partition function 1131) 
can be obtained in the large system size limit L — > oo: 



Zl{0 






1 <g< Vl + CA 
1 < Vl + CA < q. 



(18) 



one can calculate the line of canonical partition function 
zeros from 



Re g(^) =Reg 



{II) 



(24) 



in which g^^^ and g^^^^ are the free energy functions in 
the first and the second phase respectively. Using |5DJ, 
(I22|l . 123|l and l|24|l we find in the thermodynamic limit 

(L, M -^ oo, p : 



L 



For a fixed total density of particles p (which means fixed 
^) and A, the phase transition occurs at qc = y/1 + £,A. 
Now one can easily calculate the canonical partition func- 
tion of the system in these phases using l^. By using 
^ for the first phase the condition 1 < q < vT+"fA 
translates to 1 < g < . and the canonical partition 
function which is given by: 



Z 



(I) 

L.M 



2ni 



d^ 



zliO 

CM+l 



(19) 



can readily be calculated by applying the steepest decent 
method. We find 



7(1) 



(l-(l-p)<z2)(g2_(l_p)) 



,i<q< 



1 



(20) 



For the second phase the condition 1 < ^1 + ^A < q 



translates to 1 < 



/T^ 



< q. We have also 



zill) 






iM+l 



(21) 



Keeping in mind that the contour of the integral above is 
a unit circle and that its integrand has two poles, which 



one of them ^] 



_ 2_ 



is smaller than unity and the other 



^2=^1+ H is larger than unity, one can easily calculate 
(I21|l using the steepest decent method. We find 



7(11) 

'L,M 



4+2L^ 



M 



(g2 + l)(g2_l)M 



1 < 



VT 



< q- (22) 



It can be seen that for a fixed density p the transition 



pomt qc 



does not depend on A. This has al- 



ready been predicted in [l3|. For the case q < 1 the 
transition point is found to be q' = ^/l — p which agrees 
again with our predications in |l7l| . 

In [l3| we have estimated the roots of the canonical par- 
tition function Z^^m as a function q both for q > I and 
q < 1. From there we were able to find the transition 
points. Let us now calculate the line of the canonical 
partition function zeros of the model in the complex q- 
plane for q > 1. Defining the extensive part of the free 
energy as 



hm — In Zlm 

L.M^oo L 



(23) 



(1-p) 



p-i 



[{u^ 



(2ui;)2]p/2 



(25) 



in which we have defined u := Re{q) and v := Im{q). 
It can easily be verified that 1)25(1 intersects the positive 
real q-axis at u^ — ,} As can be seen the equation 

(|25|l is exactly the one that we had obtained in [2J| for 
the same model with the left boundary open and conser- 
vation of total number of particles. In [22 we had also 
found that the density of canonical partition function ze- 
ros as a function of q drops to zero near the critical point. 
This indicates that a second-order phase transition takes 
place at the critical point. We have checked that the den- 
sity of canonical partition function zeros in the present 
model also approaches to zero near the critical points qc 
and q'^. 

For q < 1 we should only change q -^ q^^ which means 
u -^ .i^ •> and V — > .7^ ■> in 1251) . In this case the line of 
canonical partition function zeros intersects the positive 
real q-axis at u'^ = ^l — p. It is not difficult to check 
that in the thermodynamic limit the numerical estimates 
for the canonical partition function zeros obtained in |l7| 
lay exactly on H25() and its counterpart for q < 1. 



IV. DENSITY PROFILE OF PARTICLES 



Now we consider the average particle number at each 
site. As for the partition functions, it turns out that 
the calculation of density profile of the particles in the 
grand canonical ensemble is much easier than that in the 
canonical ensemble; therefore, we will first calculate H10() 
and then translate out results into the canonical ensemble 
using (|11|) . The unnormalized average particle number at 
site i in the grand canonical ensemble H1(J|) can also be 
written as: 



{p^)riO^{W\C^~'^DC'^-'\V) 



(26) 



in which C :— ^D + E. Now one can use the matrix 
representation 10 to calculate (|26|) . After some algebra 
we find 



{p^)^l\0 = u,{0 q'^^-''^' + u,iO q'-'^il + CA)^- + ^3(0 q'-'\^ ^ ^^''^-' 






2L-2^^+^^^^-l 



1 



1 + A 



y-' + ueiOii + ^A) 



L-l 



+ 



(27) 



in which we have defined 



«i(0 

"2(0 
"3(0 
"4(0 
"5(0 
"6(0 
"7(0 



q\q^ - l)2eA2(e(2 + CA) - l)(g2 _ 1 _ ^A)-i 



(g2 - ^A - l)(g2(i + A) - ^A - 1)(92(1 + ^A) - l){q^l + ^A) - A - 1) ' 

-g'(g' - i)e'A 

(q2-eA-l)(g2(i+eA)-l)' 

g^(g^-l)($-l)eA 

(g2(l + A) - ^A - 1)(92(1 + ^A) - A - 1) ' 

g^g^ - l)e^A 
(g2+^A-l)(g2(i+^A)-l)' 

-gng^-i)(e-iKA 

(q2(l + A) - CA - l)(g2(i + eA) - A - 1) ' 

-(74^3^2 

(g2-eA-l)(g2(l+^A)-l)' 

g^(e-l)^eA^ 

(1 + A)(92(i + A) - CA - l)(g2(i + ^A) - A - 1) • 

I 



(28) 
(29) 
(30) 

(31) 
(32) 
(33) 
(34) 



The asymptotic behaviors of H27|) can now be distin- 
guished for the two mentioned cases. For the first case 
where 1 < g < ^/T+'^K the leading terms are the second, 
the fourth and the sixth terms in H27|) . Now using pi|) 
and H2U|I one can calculate the average particle number of 
the particles at site i in the canonical ensemble by apply- 
ing the steepest decent method. In the thermodynamic 
limit the result is: 



{p^)=p+{q^-l)[e --{l~p)e 



L-i 
«2 



I < q < 



1 



q2 



(35) 
|-i and 



in which the correlation lengths arc 1^1 = | ln(- 

^2 — I ln(q^(l — p))|^^. For a plot of this profile see figure 
2 in 17]. In the second case where 1 < ^/l 4- ^A < q 
the leading terms are the first and the fourth terms in 
(|27|l. Using numerical calculations we had predicted in 
[12] that the density profile of the particles in this phase 
is a shock in the bulk of the chain while it increases expo- 
nentially near the left boundary for g > 1. The density 
of the particles in the high-density region of the shock is 
equal to pHigh = ^ — q~^ while in the low-density region, 
it is zero plow — 0. One can easily calculate the share 



of the first term in to the density profile of the particles 
in the canonical ensemble using (|11() . By applying the 
steepest decent method one finds (1 — q~^)g^^^*. In or- 
der to calculate the share of the fourth term in (|27|l in 
the grand canonical ensemble we use the following proce- 
dure: When L is large, the average density profile can be 
described by a continuous function p(x) in terms of the 
rescaled variable x = j w here < a; < 1. By using Ijllll 
for the fourth term in l|27l) we find that the derivative of 
p{x) has the following form: 



—p{x) ~ poexp[L • F{x)\ 
ax 



where 



F{x) — —x\vLq + .T In p In ■ 



P 



X — p A(x — p) 

The constant po in (|36|l is determined by the fact that 
^1 d 



(36) 

(37) 



dx 



p{x)dx = Plow - PHigh =9 - 1- (38) 



It turns out that the function F{x) has a maximum value 



at xo 



1-9- 



One can expand F{x) around xq up to 



the second order and approximate (|36|l with a Gaussian and find: 
I 



dx 



p{x) 



L 



2'Kpq 



;il-q-^fexpi-L 



{1 - q-^)\x - XqY 
2pq-^ 



(39) 



By integrating (|39|l the average particle number at site i be 



in the canonical ensemble for 1 < 



< g is found to 



(P.) = 



= {l-q-')q'e 



^ 1 

«3 + — 



-erfc{y 



L 



2pq- 



r(l-'^-^)(T- 



L 1-q- 



r)) , 1< 



1 



VT^ 



<q 



(40) 



r 



in which the exponential part drops with the length scale 
^3 — I Ing'*!"-'^ and erfc{- • • ) is the complementary error 
function. As can be seen from H4(J|I the average particle 
number at site i far from the left boundary is an error 
function interpolating between the low-density and the 
high-density regions with width scaling as \/L. For a 
plot of this profile see figure 2 in [T3 ■ 



V. CONCLUDING REMARKS 

In this paper we studied a one-dimensional asymmetric 
branching-coalescing model with reflecting boundaries in 
a canonical ensemble where the total number of particles 
is a constant. This model has already been studied in 
literatures in a grand canonical ensemble where the total 
number of particles on the chain is not fixed and can vary 
between and 1 (see [23, Um ^^"^ references therein) . 
Without particle number conservation the parameter A, 
which is the ratio of branching to coalescing rates, gov- 
erns the average density of particles on the chain. In 
this case the phase diagram of the model consists of two 
phases: A high-density phase and a low-density phase. In 
the hight-density phase the density profile of the particles 



has an exponential behavior with two different correla- 
tion lengths |ln(Y^)|-i and \\n{q^{l + A))\-'^ . In the 
low-density phase the density profile of the particles has 
also an exponential behavior; however, with the length 
scales |ln((7^)|~^ and | ln( j-ta)|~^. On the coexistence 
line between these two phases the density profile of the 
particles has a linear decay in one end of the chain while 
it has an exponential decay in the other end of the chain 
with the length scale | ln((7^)|^-'^. 

In the canonical ensemble the total density of particles 
on the chain is controlled by the parameter p instead of 
A. With particle number conservation it turns out that 
for g > 1 the model has two different phases: An ex- 
ponential phase and a shock phase. In the exponential 
phase the density profile of the particles has an expo- 
nential behavior with two length scales | ln(:j^)|^^ and 
I \n{q^{l — p))\~^. In the shock phase the density profile of 
the particles drops exponentially near the left boundary 
with the length scale | ln(q^)|~^. In the bulk of the chain 
the density profile of the particles is an error function 
with an interface which extends over a region of width 

Vl. 
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